"""
This Script take a column file with angles and output a dictionary for each fullscan
"""
import pylab as plt
import numpy
import datetime
import cPickle
from avoscan import processing
import wx
import sys
from std_ops.iter_ import parallel_process
import scipy
import math
import os.path

#TODO list:
# Move GUI code out of this file
# create getters and setters for processed scan class
# sort out the update order of all the plots - it is a mess at the moment
# optimise redrawing to make it faster.
# dialog from choosing h_dist and v_angle and default fit order
# autofit needs to choose it's own fit order
# pull processing functions out of classes to make them easier to test/debug
# scanSelector.update_plot means that current selection is lost (not coloured red).
# skip controls should hav onnext and onprev methods and should disable buttons while loading next plot
# when one window closes, so too should the other!
#display the wangle filename in the window.

class DraggableLine:
    lock = None
    
    def __init__(self, x_pos, axis, linewidth=2.0, color='r', callback=None):
        
        self.line = axis.plot([x_pos]*2,axis.get_ylim(), linewidth=linewidth, color=color)[0]

        self.press = None
        self.background = None
        self.callback = callback
        self.connect()

   
    def connect(self):
        #connect to all the events we need
        self.cidpress = self.line.figure.canvas.mpl_connect('button_press_event', self.on_press)
        self.cidrelease = self.line.figure.canvas.mpl_connect('button_release_event', self.on_release)
        self.cidmotion = self.line.figure.canvas.mpl_connect('motion_notify_event', self.on_motion)           
        self.cidredraw = self.line.figure.canvas.mpl_connect('draw_event', self.on_redraw)

    def disconnect(self):
        #disconnect all the stored connection ids
        self.line.figure.canvas.mpl_disconnect(self.cidpress)
        self.line.figure.canvas.mpl_disconnect(self.cidrelease)
        self.line.figure.canvas.mpl_disconnect(self.cidmotion)
        self.line.figure.canvas.mpl_disconnect(self.cidredraw)
        
    def on_press(self, event):
        #on button press we will see if the mouse is over us and store some data
        if event.inaxes != self.line.axes: return
        if DraggableLine.lock is not None: return
        contains, attrd = self.line.contains(event)
        if not contains: return

        x0 = self.line.get_xdata()[0],
        self.press = x0, event.xdata, event.ydata
        DraggableLine.lock = self

        # draw everything but the selected rectangle and store the pixel buffer
        canvas = self.line.figure.canvas
        axes = self.line.axes
        self.line.set_animated(True)
        canvas.draw()
        self.background = canvas.copy_from_bbox(self.line.axes.bbox)

        # now redraw just the rectangle
        axes.draw_artist(self.line)

        # and blit just the redrawn area
        canvas.blit(axes.bbox)


    def get_xpos(self):
        return self.line.get_xdata()[0]
    
    
    def set_xpos(self,x_pos):
        self.line.set_xdata([x_pos]*2)
#        canvas = self.line.figure.canvas
#        axes = self.line.axes        
#        
#        # restore the background region
#        canvas.restore_region(self.background)
#
#        # redraw just the current rectangle
#        axes.draw_artist(self.line)
#
#        # blit just the redrawn area
#        canvas.blit(axes.bbox)
        self.line.axes.redraw_in_frame()
    
    def on_redraw(self, event):
        #reset the line's ylims so that you can't zoom to outside the size of the line
        self.line.set_ydata(self.line.axes.get_ylim())
        
        canvas = self.line.figure.canvas
        axes = self.line.axes        
        
        # restore the background region
        canvas.restore_region(self.background)

        # redraw just the current rectangle
        axes.draw_artist(self.line)

        # blit just the redrawn area
        canvas.blit(axes.bbox)
    
    
    def on_motion(self, event):
        #on motion we will move the rect if the mouse is over us
        if DraggableLine.lock is not self:
            return
        if event.inaxes != self.line.axes: return
        x0, xpress, ypress = self.press
        dx = event.xdata - xpress
        dy = event.ydata - ypress

        self.line.set_xdata([x0[0] + dx]*2)
        self.line.set_ydata(self.line.axes.get_ylim())
        self.line.set_linestyle('--')

        canvas = self.line.figure.canvas
        axes = self.line.axes
        # restore the background region
        canvas.restore_region(self.background)

        # redraw just the current rectangle
        axes.draw_artist(self.line)

        # blit just the redrawn area
        canvas.blit(axes.bbox)

 
    def on_release(self, event):
        #on release we reset the press data
        if DraggableLine.lock is not self:
            return
        
        self.line.set_linestyle('-')
        self.press = None
        DraggableLine.lock = None

        # turn off the rect animation property and reset the background
        self.line.set_animated(False)
        self.background = None

        # redraw the full figure
        self.line.figure.canvas.draw()
        
        if self.callback is not None:
            self.callback(event.xdata)


def length_of_scan_path(r, theta0, h_dist, v_angle, theta):
    """
    By assuming a cylindrical plume, this function returns the path length through the
    plume for each scan angle. For a cylindrical plume, the path through the plume
    will always be an ellipse, whose semi-minor axis is equal to the width of the plume
    (i.e. the radius of the cylinder) and whose semi-major axis is determined by the
    viewing geometry.
    
    * r - radius of the plume (of the cylinder)
    * theta0 - scan angle of the centre of the plume (in degrees)
    * h_dist - the horizontal distance of the scanner from the centre of the plume
    * v_angle - the angle (in degrees) from horizontal of the scan (e.g. if you were 
                level with the plume then this will be zero).
    * theta - numpy array of scan angles.
    
    Returns: array of path lengths through the plume (one for each scan angle). Scan angles
             that did not intersect the plume will have a path length of zero.
    
    ---------------------------------------------------------------------------------------
    A few doctests...
    
    when v_angle=0 then we are dealing with a circle
    
    >>> p = length_of_scan_path(1.0, 90.0, 10.0, 0.0, numpy.array([90.0]))
    >>> print p[0] == 2.0
    True
    >>> p = length_of_scan_path(1.0, 90.0, 10.0, 0.0, numpy.array([0.0]))
    >>> print p[0] == 0.0
    True
    >>>
    
    when v_angle=60 then we expect the semi-major axis to be 2*r
    
    >>> p = length_of_scan_path(1.0, 90.0, 10.0, 60.0, numpy.array([90.0]))
    >>> print round(p[0],4) == 4.0
    True
    >>> p = length_of_scan_path(1.0, 90.0, 10.0, 60.0, numpy.array([0.0]))
    >>> print p[0] == 0.0
    True
    >>>
    
    #when theta is such that the line is a tangent, we expect to get zero back
    #>>> p = length_of_scan_path(1.0, 90.0, 10.0, 0.0, numpy.asin
    
    ---------------------------------------------------------------------------------------
    """
    #ellipse: x^2/a^2 + y^2/b^2 = 1
    #ellipse params are constrained by the viewing angle
    a = r    
    b = r / math.cos(math.radians(v_angle))
    
    #transform h_dist into the plane of the ellipse
    h_dist = h_dist / math.cos(math.radians(v_angle))
    
    #each scan line can be defined by two points on the line - here we use the 
    #x and y axis intercepts as the points
    x1 = numpy.zeros_like(theta)
    y1 = numpy.ones_like(theta) * (-h_dist)
    
    x2 = h_dist * numpy.tan(numpy.radians(theta - theta0))
    y2 = x1.copy()
    
    #use the equation of a straight line x = cy + d we do it 
    #this way around because it is possible that theta==theta0, which 
    #leads to x1==x2 and an infinite gradient. However, h_dist is never
    #0 and so dx/dy will never be infinite.
    c = (x2 - x1) / (y2 - y1)
    d = x1 - (c * y1)
    
    #substitute the expression for x from the line into the equation of the ellipse
    #and solve using the quadratic equation.
    #quadratic equation -> x = (-B +- sqrt(B^2 - 4AC))/(2A)
    A = (a ** 2 / b ** 2) + c ** 2
    B = 2.0 * c * d
    C = d ** 2 - a ** 2

    #mask out scan angles that did not intersect the plume (they have imaginary roots)        
    discriminant = B ** 2 - 4 * A * C
    invalid_values = discriminant <= 0
    
    A = numpy.ma.masked_where(invalid_values, A, copy=False)
    B = numpy.ma.masked_where(invalid_values, B, copy=False)
    C = numpy.ma.masked_where(invalid_values, C, copy=False)
    c = numpy.ma.masked_where(invalid_values, c, copy=False)
    d = numpy.ma.masked_where(invalid_values, d, copy=False)

    #solve for the intersection points of the lines with the ellipse    
    cy1 = (-B + numpy.sqrt(B ** 2 - 4.0 * A * C)) / (2.0 * A)
    cy2 = (-B - numpy.sqrt(B ** 2 - 4.0 * A * C)) / (2.0 * A)
        
    cx1 = c * cy1 + d
    cx2 = c * cy2 + d
    
    #return the distances between the intersection points (i.e. the path lengths
    #and fill in the values that did not intersect with zeros.
    return numpy.ma.filled(numpy.ma.sqrt((cx2 - cx1) ** 2 + (cy2 - cy1) ** 2), 0.0)
    


def estimate_plume_radius(angles, ydata, h_dist, v_angle, subplot=None):
    """
    Assuming a cylindrical plume, this function finds the plume radius and 
    position (scan angle to the centre) that best fits the supplied data.
    It works by doing a least squares fit of linearly scaled path lengths
    through the plume for different radii to the data. In this way it is 
    assuming that the plume is uniform and that changes in column amounts
    are only due to different path lengths.
    
    * angles - numpy array of scan angles
    * ydata - numpy array of column amounts or ratios or similar
    * h_dist - the horizontal distance of the scanner from the centre of the plume
    * v_angle - the angle (in degrees) from horizontal of the scan (e.g. if you were 
                level with the plume then this will be zero).
    * subplot - a matplotlib Axes object that you want the best fit function to be
                plotted into, or None if you don't want anything plotted.
                
    Returns: a tuple of best fit parameters (radius, theta0) where theta0 is the scan
             angle towards the centre of the plume.
    
    ---------------------------------------------------------------------------------------
    A few doctests...
    
    >>> ang = numpy.array(range(45, 136))
    >>> so2 = length_of_scan_path(200.0, 90.0, 1000.0, 45, ang)
    >>> r,theta0 = estimate_plume_radius(ang, so2, 1000.0, 45)
    >>> print round(r,3) == 200.0
    True
    >>> print round(theta0,3) == 90.0
    True
    
    """

    # guess some fit parameters - we use the width at half height as an estimate of radius
    ymax = ydata.max()

    #print "weighted av = ", numpy.average(angles,weights=numpy.abs(ydata))

    angle_at_half_height = angles[numpy.abs(ydata - (ymax / 2.0)).argmin()]
       
    theta0_guess = numpy.average(angles,weights=numpy.abs(ydata))
    
    variance = numpy.dot(numpy.abs(ydata), (angles-theta0_guess)**2)/numpy.abs(ydata).sum()  # Fast and numerically precise
    
    r_guess_angle = math.sqrt(variance)
    r_guess = h_dist * math.tan(math.radians(r_guess_angle))
    
    #print "theta0 guess = ", theta0_guess
    
    #r_guess = h_dist * math.tan(math.radians(math.fabs(theta0_guess - angle_at_half_height)))
    #print "r_guess = ",r_guess
    scaling_guess = ymax / (2.0 * (r_guess / math.cos(math.radians(v_angle))))

    #put guess params into an array ready for fitting
    p0 = numpy.array([r_guess, theta0_guess, scaling_guess])

    #define the path length function and associated error function
    fitfunc = lambda p, x: p[2] * length_of_scan_path(p[0], p[1], h_dist, v_angle, x)
    errfunc = lambda p, x, y: fitfunc(p, x) - y
   
    # do the fitting
    p1, success = scipy.optimize.leastsq(errfunc, p0, args=(angles, ydata))

    if success not in (1, 2, 3, 4):
        raise RuntimeError, "Could not fit cylindrical plume to data."

    if subplot is not None:
        subplot.plot(angles, ydata, 'g+', label='data')
        
        #plot the ellipse fit at higher angle resolution
        high_res_ang = numpy.linspace(numpy.min(angles), numpy.max(angles), 1000)
        subplot.plot(high_res_ang, fitfunc(p1, high_res_ang) , 'r-', label='fit')
    
    return (p1[0], p1[1])



def find_residuals(x, y, x_err, y_err, func):
    """
    Calculates the absolute distances of a data series from a trend line.
    However, any distances that are within the error bars are set 
    to zero. This makes it easy to see if a trend line fits some data within
    experimental error. 
    
    * x - numpy array of x data
    * y - numpy array of corresponding y data
    * x_err - absolute error in x direction (x +- x_err)
    * y_err - absolute error in y direction (y +- y_err)
    * func - a Python function that describes the trendline such that y_trend = func(x)
             where x is a single x value.
    """
    #convert func to numpy ufunc
    fit_func = numpy.frompyfunc(func, 1, 1)
    
    x_mins = x - x_err
    x_maxs = x + x_err
    
    #evaluate the distance between the data and the trendline at the extremes of
    #the x error.
    resid_mins = numpy.abs(y - fit_func(x_mins))
    resid_maxs = numpy.abs(y - fit_func(x_maxs))
    
    #get the smallest distance of each point from the fit line
    resid = numpy.minimum(resid_mins, resid_maxs)
    
    #set residuals that are within the error bars to zero
    resid[numpy.where(resid < y_err)] = 0
    
    return resid
    


def calculate_bkgd(angles, col_amount, angle_err, col_err, fit_order=2, subplot_limits=None, subplot_calc=None):
    """
    Calculates a polynomial fit to the background for a single scan. This is done
    by fitting a Gaussian to the scan, and treating all points 2.5*sigma away from the
    mean as background points. Points which cannot be fitted within error are then removed
    one by one and the fitting is repeated each time.
    
    * angles - numpy array of scan angles
    * col_amount - numpy array of column amounts or ratios or similar
    * angle_err - numpy array of absolute errors in scan angles
    * col_err - numpy array of absolute errors in the column amounts
    * fit_order - order of polynomial that you want to fit
    * subplot - a matplotlib Axes object that you want the background fit function to be
                plotted into, or None if you don't want anything plotted.
    
    Returns: a numpy array of background levels (with the same length as the supplied 
             col_amount array).
    """
    
    #Do the Gaussian fitting
    gfit_params = processing.fit_gaussian(angles, col_amount)
       
    #mask the angles and col_amounts array to only expose the out of plume sections
    a = 2.5
    is_lower_bkgd = angles < (gfit_params.mean - a * gfit_params.sigma)
    is_higher_bkgd = angles > (gfit_params.mean + a * gfit_params.sigma)
    bkgd_mask = numpy.logical_or(is_lower_bkgd, is_higher_bkgd)
    bkgd_idxs = numpy.where(bkgd_mask)
    
    #fit a polynomial through the background
    ignored_pts = []
    
    while True:
        bkgd_angles = angles[bkgd_idxs]
        bkgd_col_amounts = col_amount[bkgd_idxs]
        bkgd_angle_err = angle_err[bkgd_idxs]
        bkgd_col_err = col_err[bkgd_idxs]
        
        polyfit_params = numpy.polyfit(bkgd_angles, bkgd_col_amounts, fit_order)
        bkgd_function = numpy.poly1d(polyfit_params)
             
        #calculate the residuals - ignoring those that are within error
        resid = find_residuals(bkgd_angles, bkgd_col_amounts, bkgd_angle_err, bkgd_col_err, bkgd_function)
        
        if len(numpy.nonzero(resid)[0]) == 0:
            break
        
        #remove the point with the greatest residual from the fit
        idx_to_skip = bkgd_idxs[0][numpy.argmax(resid)]
        ignored_pts.append((angles[idx_to_skip], col_amount[idx_to_skip], angle_err[idx_to_skip], col_err[idx_to_skip]))
        bkgd_idxs = (bkgd_idxs[0][numpy.where(bkgd_idxs[0] != idx_to_skip)],)
    
    print "Ignored", len(ignored_pts), "points from background fitting."
   
    #if we have been passed a subplot object - then plot the fit into it     
    if subplot_calc is not None:
        subplot_limits.errorbar(angles, col_amount, col_err, angle_err, '+')
        globals()["d"] = DraggableLine(gfit_params.mean - a * gfit_params.sigma, subplot_limits)
        globals()["d2"] = DraggableLine(gfit_params.mean + a * gfit_params.sigma, subplot_limits)
        #subplot_limits.plot([gfit_params.mean - a * gfit_params.sigma]*2,subplot_limits.get_ylim(),'b-', linewidth=2)
        #subplot_limits.plot([gfit_params.mean + a * gfit_params.sigma]*2,subplot_limits.get_ylim(),'g-', linewidth=2)
        subplot_limits.grid()
         
        not_bkgd_idxs = numpy.where(numpy.logical_not(bkgd_mask))
        not_bkgd_angles = angles[not_bkgd_idxs]
        not_bkgd_col_amounts = col_amount[not_bkgd_idxs]
        not_bkgd_col_err = col_err[not_bkgd_idxs]
        not_bkgd_angle_err = angle_err[not_bkgd_idxs]
        
        ignored = zip(*ignored_pts)
        if ignored:
            subplot_calc.errorbar(ignored[0], ignored[1], ignored[3], ignored[2], 'r+')
        
        #subplot_calc.errorbar(not_bkgd_angles,not_bkgd_col_amounts, not_bkgd_col_err, not_bkgd_angle_err,'g+')        
        subplot_calc.errorbar(bkgd_angles, bkgd_col_amounts, bkgd_col_err, bkgd_angle_err, 'b+')
        subplot_calc.plot(angles, [bkgd_function(i) for i in angles], 'm-', linewidth=2)
    
    #return the calculated background for all angles
    return numpy.array([bkgd_function(i) for i in angles])

def integrate(angles, data, radius, h_dist, v_angle, theta0):
    """
    
    
    ---------------------------------------------------------------------------------------
    A few doctests...
    
    >>> ang = numpy.array(range(45, 136))
    >>> so2 = length_of_scan_path(200.0, 90.0, 1000.0, 45, ang)
    >>> i = integrate(ang, so2, 200.0, 1000.0, 45.0, 90.0)
    >>> print round(i, 4) == round(200**2 * math.pi, 4)
    True
  
    """
    #TODO - data needs to be converted to something per m^2 rather than whatever it is in currently
    #TODO - error analysis!
    
    #the integration assumes the angles are in increasing order so sort them here
    c = numpy.rec.fromarrays([angles, data])
    c.sort() #sort based on angles
    angles = c.f0
    data = c.f1

        
    #calculate the path length through the plume for each of the measurements
    measured_dists = length_of_scan_path(radius, theta0, h_dist, v_angle, angles)
    
    #crop data to only be working inside the plume
    pts_inside_plume = numpy.where(measured_dists > 0)
    angles = angles[pts_inside_plume]
    data = data[pts_inside_plume]
    measured_dists = measured_dists[pts_inside_plume]
    
    #calculate the areas of the circle sections between each scan line
    #for chord area calculation see: http://mathworld.wolfram.com/Chord.html
    r = h_dist * numpy.sin(numpy.radians(angles - theta0))    
    chord_areas = numpy.zeros((len(angles) + 1,), dtype="float")
    chord_areas[-1] = (math.pi * radius ** 2)
    chord_areas[:-1] = chord_areas[-1] - ((radius ** 2 * numpy.arccos(r / radius)) - (r * numpy.sqrt(radius ** 2 - r ** 2)))
    
    #scan areas are equal to the difference between chord areas
    scan_areas = numpy.zeros_like(chord_areas)
    scan_areas[0] = chord_areas[0]
    scan_areas[1:] = chord_areas[1:] - chord_areas[:-1]
       
    #calculate the averaged concentraions
    concs = data / measured_dists
    
    av_concs = numpy.zeros_like(scan_areas)
    av_concs[0] = concs[0] #don't bother averaging at the edges
    av_concs[-1] = concs[-1] #don't bother averaging at the edges
    av_concs[1:-1] = (concs[:-1] + concs[1:]) / 2.0 
    
    
    #calculate the integral
    return numpy.sum(av_concs * scan_areas)
    

class ProcesssedScan:
    def __init__(self, angles, times, col_amounts, angle_err, col_err, h_dist, v_angle, bkgd_fit_order=1, bkgd_threshold=5.0):
        self.__is_calculated = False
        
        self.angles = angles
        self.times = times    
        self.col_amount = col_amounts
        self.angle_err = angle_err
        self.col_err = col_err
        self.h_dist = h_dist
        self.v_angle = v_angle
        
        #background fitting parameters
        self.bkgd_fit_order = bkgd_fit_order
        self.bkgd_threshold = bkgd_threshold
        self.manual_bkgd = False
    
    def get_bkgd_subtracted_plume_points(self):
        #TODO - what if this stuff hasn't been calculated yet?
        plume_mask = numpy.logical_not(self.bkgd_mask)
        plume_idxs = numpy.where(plume_mask)        
        bkgd_subtracted_ydata = self.col_amount - self.bkgd_func(self.angles)
        return self.angles[plume_idxs], self.col_amount[plume_idxs], self.angle_err[plume_idxs], self.col_err[plume_idxs]
        
                
    def calculate(self):
        print self.times[0]
        self.__calculate_bkgd()
        self.__estimate_ica()
        self.__is_calculated = True
        return self
    
    def plot_radius_fit(self, ax):
        plume_mask = numpy.logical_not(self.bkgd_mask)
        plume_idxs = numpy.where(plume_mask)
        if self.ica == 0:
            ax.clear()
            return
        bkgd_subtracted_ydata = self.col_amount[plume_idxs] - self.bkgd_func(self.angles)[plume_idxs]
        ax.plot(self.angles[plume_idxs]+self.angle_err[plume_idxs], bkgd_subtracted_ydata + self.col_err[plume_idxs], 'g-', linewidth=0.5)
        ax.plot(self.angles[plume_idxs]+self.angle_err[plume_idxs], bkgd_subtracted_ydata- self.col_err[plume_idxs], 'g-', linewidth=0.5)
        ax.plot(self.angles[plume_idxs]-self.angle_err[plume_idxs], bkgd_subtracted_ydata + self.col_err[plume_idxs], 'g-', linewidth=0.5)
        ax.plot(self.angles[plume_idxs]-self.angle_err[plume_idxs], bkgd_subtracted_ydata - self.col_err[plume_idxs], 'g-', linewidth=0.5)
        
        #plot the ellipse fit at higher angle resolution
        high_res_ang = numpy.linspace(numpy.min(self.angles[plume_idxs]-self.angle_err[plume_idxs]), numpy.max(self.angles[plume_idxs]+self.angle_err[plume_idxs]), 1000)
        
        for radius, theta0 in self.plume_radii:
            ax.plot(high_res_ang,length_of_scan_path(radius, theta0, self.h_dist, self.v_angle, high_res_ang), 'r-')
    
    
    def plot_bkgd_fit(self, ax):
        
        ignored = zip(*self.bkgd_ignored_pts)
        if ignored:
            ax.errorbar(ignored[0], ignored[1], ignored[3], ignored[2], 'r+')
        
        #subplot_calc.errorbar(not_bkgd_angles,not_bkgd_col_amounts, not_bkgd_col_err, not_bkgd_angle_err,'g+')        
        ax.errorbar(self.angles[self.bkgd_idxs], self.col_amount[self.bkgd_idxs], self.col_err[self.bkgd_idxs], self.angle_err[self.bkgd_idxs], 'b+')
        ax.plot(self.angles, [self.bkgd_func(i) for i in self.angles], 'm-', linewidth=2)
   
   
      
    def __calculate_bkgd(self):
        
        if not self.manual_bkgd:
            #Do the Gaussian fitting in order to estimate the width of the plume and it's position
            if len(self.angles) != len(self.col_amount):
                print "lengths do not match"
            gfit_params = processing.fit_gaussian(self.angles, self.col_amount)
           
            #mask the angles and col_amounts array to only expose the out of plume sections
            self.bkgd_lowlim = gfit_params.mean - math.fabs(self.bkgd_threshold * gfit_params.sigma)
            self.bkgd_highlim = gfit_params.mean + math.fabs(self.bkgd_threshold * gfit_params.sigma)

            if self.bkgd_lowlim > self.bkgd_highlim or self.bkgd_lowlim < min(self.angles[0],self.angles[-1]): #remember scans can be in either direction
                #this quite often happens if there are insufficient points in the plume to effectively
                #fit a gaussian to it. As a hack - we just reset the lims to include all of the data
                # there is probably a smarter way to do this
                print "time = ", self.times[0], "mean = ", gfit_params.mean, "sigma = ", gfit_params.sigma
                self.bkgd_lowlim = gfit_params.mean
                self.bkgd_highlim = gfit_params.mean
        
        is_lower_bkgd = self.angles < self.bkgd_lowlim
        is_higher_bkgd = self.angles > self.bkgd_highlim
        self.bkgd_mask = numpy.logical_or(is_lower_bkgd, is_higher_bkgd)
        bkgd_idxs = numpy.where(self.bkgd_mask)
        
        #fit a polynomial through the background
        ignored_pts = []
        
        while True:
            bkgd_angles = self.angles[bkgd_idxs]
            bkgd_col_amounts = self.col_amount[bkgd_idxs]
            bkgd_angle_err = self.angle_err[bkgd_idxs]
            bkgd_col_err = self.col_err[bkgd_idxs]
            
            polyfit_params = numpy.polyfit(bkgd_angles, bkgd_col_amounts, self.bkgd_fit_order)
            bkgd_function = numpy.poly1d(polyfit_params)
                 
            #calculate the residuals - ignoring those that are within error
            resid = find_residuals(bkgd_angles, bkgd_col_amounts, bkgd_angle_err, bkgd_col_err, bkgd_function)
            
            if len(numpy.nonzero(resid)[0]) == 0:
                break
            
            #remove the point with the greatest residual from the fit
            idx_to_skip = bkgd_idxs[0][numpy.argmax(resid)]
            ignored_pts.append((self.angles[idx_to_skip], self.col_amount[idx_to_skip], self.angle_err[idx_to_skip], self.col_err[idx_to_skip]))
            bkgd_idxs = (bkgd_idxs[0][numpy.where(bkgd_idxs[0] != idx_to_skip)],)
        
        #store the calculated background for all angles
        self.bkgd_func = bkgd_function
        self.bkgd_ignored_pts = ignored_pts
        self.bkgd_idxs = bkgd_idxs

        
    def __estimate_ica(self):
        
        plume_mask = numpy.logical_not(self.bkgd_mask)
        plume_idxs = numpy.where(plume_mask)
        
        bkgd_subtracted_ydata = self.col_amount - self.bkgd_func(self.angles)
        
        #do the calculation for the full range of errors
        tmp = [(self.angles[plume_idxs] + self.angle_err[plume_idxs], bkgd_subtracted_ydata[plume_idxs] + self.col_err[plume_idxs]),
               (self.angles[plume_idxs] + self.angle_err[plume_idxs], bkgd_subtracted_ydata[plume_idxs] - self.col_err[plume_idxs]),
               (self.angles[plume_idxs] - self.angle_err[plume_idxs], bkgd_subtracted_ydata[plume_idxs] + self.col_err[plume_idxs]),
               (self.angles[plume_idxs] - self.angle_err[plume_idxs], bkgd_subtracted_ydata[plume_idxs] - self.col_err[plume_idxs])]

        tmp_ica  =[]
        measurement_times=[]
        time_errors=[]
        self.plume_radii = []
        for ang, so2 in tmp:
            if len(so2) < 5:
                self.ica = 0
                self.ica_err = 0
                min_time = min(self.times)
                max_time = max(self.times)
                self.ica_time = min_time + ((min_time - max_time) / 2)
                self.ica_time_err = abs(self.ica_time - min_time)
                return

            #first of all we estimate the plume radius
            # guess some fit parameters - we use the width at half height as an estimate of radius
            ymax = so2.max()
        
            angle_at_half_height = ang[numpy.abs(so2 - (ymax / 2.0)).argmin()]
               
            theta0_guess = numpy.average(ang,weights=numpy.abs(so2))
            
            variance = numpy.dot(numpy.abs(so2), (ang-theta0_guess)**2)/numpy.abs(so2).sum()  # Fast and numerically precise
            
            r_guess_angle = math.sqrt(variance)
            r_guess = self.h_dist * math.tan(math.radians(r_guess_angle))
    
            scaling_guess = ymax / (2.0 * (r_guess / math.cos(math.radians(self.v_angle))))
        
            #put guess params into an array ready for fitting
            p0 = numpy.array([r_guess, theta0_guess, scaling_guess])
        
            #define the path length function and associated error function
            fitfunc = lambda p, x: p[2] * length_of_scan_path(p[0], p[1], self.h_dist, self.v_angle, x)
            errfunc = lambda p, x, y: fitfunc(p, x) - y
           
            # do the fitting
            p1, success = scipy.optimize.leastsq(errfunc, p0, args=(ang, so2))
        
            if success not in (1, 2, 3, 4):
                raise RuntimeError, "Could not fit cylindrical plume to data."

        
            self.plume_radii.append((p1[0], p1[1]))
            tmp_ica.append(integrate(ang, so2, p1[0], self.h_dist, self.v_angle, p1[1]))
            
            idx_theta0 = numpy.abs(ang - p1[1]).argmin()
            measurement_times.append(self.times[idx_theta0])
            
            entry_angle = p1[1] - math.degrees(math.atan(p1[0] / self.h_dist))
            idx_entry_angle = numpy.abs(ang - entry_angle).argmin()
        
            time_errors.append(abs(self.times[idx_theta0] - self.times[idx_entry_angle]))
            
            
        #find min and max ica with error
        min_ica = min(tmp_ica)
        max_ica = max(tmp_ica)
        self.ica = (min_ica + max_ica) / 2.0
        self.ica_err = abs(self.ica - min_ica)
        
        min_time = min(measurement_times)
        max_time = max(measurement_times)
        self.ica_time = min_time + ((min_time - max_time) / 2)
        self.ica_time_err = abs(self.ica_time - min_time)


def estimate_ica_error(scan, n_iter):
    
    angles, col_amount, angle_err, col_err = scan.get_bkgd_subtracted_plume_points()
    
    #create arrays of random errors
    rnd_nums_angles = numpy.random.random_sample(size=(n_iter,len(angles)))
    rnd_nums_angles = rnd_nums_angles*(2*angle_err)-angle_err
    rnd_nums_col_amounts = numpy.random.random_sample(size=(n_iter,len(col_amount)))
    rnd_nums_col_amounts = rnd_nums_col_amounts*(2*col_err)-col_err
    
    
    #add the random errors to our data
    angles = rnd_nums_angles + angles
    col_amount = rnd_nums_col_amounts + col_amount
      
    #make a first guess at fit parameters
    ymax = col_amount.max(axis=1)

    theta0_guess = numpy.average(angles, axis=1, weights=numpy.abs(col_amount))
    
    #TODO - do variance calculation as an array operation
    variance = numpy.zeros((n_iter,1), dtype=float)
    for i in range(n_iter):                         
        variance[i] = numpy.dot(numpy.abs(col_amount[i]), (angles[i]-theta0_guess[i])**2)/numpy.abs(col_amount[i]).sum()  # Fast and numerically precise
            
    r_guess_angle = numpy.sqrt(variance)
    r_guess = (scan.h_dist * numpy.tan(numpy.radians(r_guess_angle))).flatten()


    scaling_guess = ymax / (2.0 * (r_guess / math.cos(math.radians(scan.v_angle))))

    #put guess params into an array ready for fitting
    p0 = numpy.array([r_guess, theta0_guess, scaling_guess]).swapaxes(1,0)
    
    #define the path length function and associated error function
    fitfunc = lambda p, x: p[2] * length_of_scan_path(p[0], p[1], scan.h_dist, scan.v_angle, x)
    errfunc = lambda p, x, y: fitfunc(p, x) - y     
    
    # do the fitting
    p1 = numpy.zeros_like(p0)
    globals()["myi"] = 0 
    res = numpy.apply_along_axis(myfunc, 1, p0,scan,angles, col_amount)
    print res.shape
#    success = numpy.zeros((n_iter,1), dtype='float')
#    for i in range(n_iter):
#        p1[i], success[i] = scipy.optimize.leastsq(errorfunc, p0[i], args=(scan,angles, col_amount))
#        globals()["myi"]+=1
#    print "Average plume radius = ",numpy.mean(p1[:,0].T)
#    print "Std dev = ",numpy.std(p1[:,0].T)
import cProfile     
def errorfunc(p,scan,x,y):
    return (p[2] * length_of_scan_path(p[0], p[1], scan.h_dist, scan.v_angle, x[globals()["myi"]]))-y[globals()["myi"]]

def myfunc(p,scan,x,y):
    r = scipy.optimize.leastsq(errorfunc,p,args=(scan,x,y))
    globals()["myi"]+=1
    return r
    
class ScanSelector:
    def __init__(self, scans_list, ax, callback=None):
        self.ax = ax
        self.plotted_pts = {}
        self.callback = callback
        self.scans = {}
        for s in scans_list:
            self.scans[s.ica_time] = s
        
        self.times = numpy.array([s.ica_time for s in scans_list])
        self.times.sort()
        self.previous_selection = None
        self.update_plot()
        #start off with the first point selected
        self.select(self.scans[numpy.min(self.times)])

        ax.figure.canvas.mpl_connect('pick_event', self.onpick)
        #ax.set_title("Integrated Column Amount")
    
    
    def update_key_for_scan(self, scan, old_key):
        self.scans.pop(old_key)
        self.scans[scan.ica_time] = scan
        self.times[numpy.where(self.times == old_key)] = scan.ica_time
        v = self.plotted_pts.pop(old_key)
        self.plotted_pts[scan.ica_time] = v
    
    def update_plot(self):
        self.ax.clear()
        self.plotted_pts = {}
        for s in self.scans.values():
            self.plotted_pts[s.ica_time] = self.ax.errorbar(s.ica_time, s.ica, s.ica_err, s.ica_time_err, 'b.',picker=3)
        if self.previous_selection is not None:
            self.select(self.previous_selection, force_redraw=True)
        
    def delete_current_scan(self):
        self.scans.pop(self.previous_selection.ica_time)
        self.times = self.times[numpy.where(self.times != self.previous_selection.ica_time)]
        self.previous_selection = None
        self.update_plot()

        self.select(self.scans[numpy.min(self.times)])
        
    def get_current_selection(self):
        return self.previous_selection
    
    def get_current_scan_num(self):
        #returns a tuple (cur_scan_num, n) where n is the total number of scans
        n = len(self.times)
        
        cur_scan_num = numpy.where(self.times == self.previous_selection.ica_time)[0][0] + 1 #want to start counting at 1
        return (cur_scan_num, n)
        
    def onpick(self, event):
        x_pos = event.artist.get_xdata()[0]
        scan = self.scans[x_pos]
        self.select(scan)
        
    
    def select_next(self):
        current_time_idx = numpy.where(self.times == self.previous_selection.ica_time)[0][0]
        if current_time_idx + 1 == len(self.times):
            #go back to the beginning
            next_scan = self.scans[self.times[0]]
        else:
            next_scan = self.scans[self.times[current_time_idx + 1]]
        self.select(next_scan)
    
    def select_prev(self):
        current_time_idx = numpy.where(self.times == self.previous_selection.ica_time)[0][0]
        if current_time_idx == 0:
            #go to the end
            prev_scan = self.scans[self.times[-1]]
        else:
            prev_scan = self.scans[self.times[current_time_idx - 1]]
        self.select(prev_scan)
                
    def select(self, scan, force_redraw=False):
        if scan == self.previous_selection and not force_redraw:
            return
        
        plt_pt = self.plotted_pts[scan.ica_time]
        plt_pt[0].set_color('r')
        for p in plt_pt[1]:
            p.set_color('r')
        for p in plt_pt[2]:
            p.set_color(((1,0,0),))
               
        #change old selection back to blue
        if self.previous_selection is not None and scan != self.previous_selection:
            old_pt = self.plotted_pts[self.previous_selection.ica_time]
            old_pt[0].set_color('b')
            for p in old_pt[1]:
                p.set_color('b')
            for p in old_pt[2]:
                p.set_color(((0,0,1),))
        
        plt_pt[0].figure.canvas.draw()
        self.previous_selection = scan
        
        #run the callback function
        if self.callback is not None:
            self.callback(scan)


class BackgroundRangeSelector:
    def __init__(self, ax):
        self.ax = ax
        self.lim_line1 = None
        self.lim_line2 = None
        self.plot = []
    
    def update(self, scan):
        if self.plot:
            self.plot[0].remove()
            for p in self.plot[1]:
                p.remove()
            for p in self.plot[2]:
                p.remove()
        if self.lim_line1 is None:
                self.lim_line1 = DraggableLine(scan.bkgd_lowlim, self.ax, color='g')
                self.lim_line2 = DraggableLine(scan.bkgd_highlim, self.ax, color='g')
        else:
            self.lim_line1.set_xpos(scan.bkgd_lowlim)
            self.lim_line2.set_xpos(scan.bkgd_highlim)

        self.plot = self.ax.errorbar(scan.angles, scan.col_amount, scan.col_err, scan.angle_err,'b+')
        
    def get_lims(self):
        pos1 = self.lim_line1.get_xpos()
        pos2 = self.lim_line2.get_xpos()
        
        return min(pos1,pos2), max(pos1,pos2)

class BackgroundFitDisplay:
    def __init__(self, ax):
        self.ax = ax
    
    def update(self, scan):
        self.ax.clear()
        scan.plot_bkgd_fit(self.ax)

class RadiusFitDisplay:
    def __init__(self, ax):
        self.ax = ax
    
    def update(self, scan):
        self.ax.clear()
        scan.plot_radius_fit(self.ax)        

class PlotManager:
    def __init__(self,scans):
        self.callbacks = []
        #create the plotting window
        bkgd_select_ax = plt.subplot2grid((3,2), (0,0), colspan=2)
        bkgd_fit_ax = plt.subplot2grid((3,2), (1,0))
        rad_est_ax = plt.subplot2grid((3,2), (1, 1))
        ica_ax = plt.subplot2grid((3,2), (2, 0), colspan=2)
        
        self.bkgd_selector = BackgroundRangeSelector(bkgd_select_ax)
        self.bkgd_fit = BackgroundFitDisplay(bkgd_fit_ax)
        self.rad_fit = RadiusFitDisplay(rad_est_ax)
        self.scan_selector = ScanSelector(scans, ica_ax, callback=self.update)
    
    
    def delete_current_scan(self):
        self.scan_selector.delete_current_scan()    
        
    def register_callback(self, f):
        self.callbacks.append(f)
        
    def get_scans(self):
        return self.scan_selector.scans.values()
    
    def next(self, scan):
        self.scan_selector.select_next()
        
    def prev(self, scan):
        self.scan_selector.select_prev()
    
    def show(self):    
        plt.show()
    
    def get_current_scan_num(self):
        #returns a tuple (cur_selection, n) where n is the total number of scans
        return self.scan_selector.get_current_scan_num() 
    
    def get_current_scan(self):
        return self.scan_selector.get_current_selection()
             
    def update(self, scan, old_key=None):
        #if the updated scan has a new ica_time, then we need to update the dict
        if old_key is not None:
            self.scan_selector.update_key_for_scan(scan, old_key)
        
        wx.BeginBusyCursor()
        self.bkgd_selector.update(scan)
        self.bkgd_fit.update(scan)
        self.rad_fit.update(scan)
        
        plt.draw()
        
        for f in self.callbacks:
            f(scan)

        wx.EndBusyCursor()


class ControlWindow(wx.Frame):
    def __init__(self, plot_manager, title):
        self.plot_manager = plot_manager
        wx.Frame.__init__(self, None, wx.ID_ANY, title)
        self.top_panel = wx.Panel(self, wx.ID_ANY)
        self.main_sizer = wx.BoxSizer(wx.VERTICAL)
        
        self.main_sizer.Add(SkipControls(self.top_panel, self.plot_manager),0,wx.ALIGN_TOP|wx.ALIGN_CENTER_HORIZONTAL)
        self.main_sizer.Add(BackgroundFittingControls(self.top_panel, self.plot_manager),0,wx.EXPAND)
        
        self.hsizer = wx.BoxSizer(wx.HORIZONTAL)
        
        self.save_all_button = wx.Button(self.top_panel, wx.ID_ANY, "Save All")
        self.delete_button = wx.Button(self.top_panel, wx.ID_ANY, "Delete Scan")
        self.hsizer.Add(self.delete_button,0,wx.ALIGN_LEFT|wx.ALIGN_BOTTOM)
        self.hsizer.Add(self.save_all_button,0,wx.ALIGN_RIGHT|wx.ALIGN_BOTTOM)
        self.main_sizer.Add(self.hsizer,1,wx.ALIGN_BOTTOM|wx.ALIGN_RIGHT)
        
        self.top_panel.SetSizer(self.main_sizer)
        self.main_sizer.Fit(self)
        self.top_panel.SetAutoLayout(1)
        
        wx.EVT_BUTTON(self, self.save_all_button.GetId(), self.on_save)
        wx.EVT_BUTTON(self, self.delete_button.GetId(), self.on_delete)
        self.Show()
    
    def on_save(self,evt):
        filename = wx.SaveFileSelector("Save","scans")
        wx.BeginBusyCursor()
        if not filename.endswith(".scans"):
            filename = filename + ".scans"
        with open(filename, "wb") as ofp:
            cPickle.dump(self.plot_manager.get_scans(), ofp)
        wx.EndBusyCursor()
            
    def on_delete(self,evt):
        self.plot_manager.delete_current_scan()
        

class  SkipControls(wx.Panel):
    def __init__(self, parent, plot_manager):
        wx.Panel.__init__(self, parent, wx.ID_ANY)
        self.plot_manager = plot_manager
        self.vsizer = wx.BoxSizer(wx.VERTICAL)
        self.hsizer = wx.BoxSizer(wx.HORIZONTAL)
        
        s,n = self.plot_manager.get_current_scan_num()
        self.txt = wx.StaticText(self,wx.ID_ANY,"Displaying scan "+str(s)+" of "+str(n))
        
        self.vsizer.AddSpacer(10)
        self.vsizer.Add(self.txt,0,wx.ALIGN_CENTER_HORIZONTAL)
        
        self.next_button = wx.Button(self, wx.ID_ANY, ">")
        self.prev_button = wx.Button(self, wx.ID_ANY, "<")
        
        self.hsizer.Add(self.prev_button,0,wx.ALIGN_LEFT)
        self.hsizer.Add(self.next_button,0,wx.ALIGN_RIGHT)
        
        self.vsizer.AddSpacer(5)        
        self.vsizer.Add(self.hsizer, 1, wx.ALIGN_CENTER_HORIZONTAL|wx.EXPAND)
        self.vsizer.AddSpacer(10)    
        wx.EVT_BUTTON(self, self.next_button.GetId(), self.plot_manager.next)
        wx.EVT_BUTTON(self, self.prev_button.GetId(), self.plot_manager.prev)
        
        self.plot_manager.register_callback(self.update)
        
        self.SetSizer(self.vsizer)
        self.vsizer.Fit(self)
        self.SetAutoLayout(1)
            
    def update(self, scan):
        s,n = self.plot_manager.get_current_scan_num()
        self.txt.SetLabel("Displaying scan "+str(s)+" of "+str(n))
        

class BackgroundFittingControls(wx.Panel):
    def __init__(self, parent, plot_manager):
        self.plot_manager = plot_manager
        wx.Panel.__init__(self, parent, wx.ID_ANY, style=wx.BORDER_SUNKEN)
        
        self.vsizer = wx.BoxSizer(wx.VERTICAL)
        self.autofit_check_box = wx.CheckBox(self, wx.ID_ANY, "Auto fit")
        
        self.vsizer.Add(self.autofit_check_box,0,wx.ALIGN_LEFT|wx.ALIGN_TOP)
        
        
        self.lims_hsizer = wx.BoxSizer(wx.HORIZONTAL)
        self.lower_lim_txtbox = wx.TextCtrl(self, wx.ID_ANY)
        self.upper_lim_txtbox = wx.TextCtrl(self, wx.ID_ANY)
        self.get_lims_button = wx.Button(self, wx.ID_ANY, "Get Limits")
        
        
        self.lims_hsizer.AddSpacer(5)
        self.lims_hsizer.Add(wx.StaticText(self, wx.ID_ANY,"Lower Limit:"),0, wx.ALIGN_CENTER_VERTICAL)
        self.lims_hsizer.Add(self.lower_lim_txtbox,0, wx.ALIGN_CENTER_VERTICAL)
        self.lims_hsizer.AddSpacer(10)
        self.lims_hsizer.Add(wx.StaticText(self, wx.ID_ANY,"Upper Limit:"),0, wx.ALIGN_CENTER_VERTICAL)
        self.lims_hsizer.Add(self.upper_lim_txtbox,0, wx.ALIGN_CENTER_VERTICAL)
        self.lims_hsizer.AddSpacer(10)
        self.lims_hsizer.Add(self.get_lims_button,0, wx.ALIGN_CENTER_VERTICAL)
        self.lims_hsizer.AddSpacer(10)
        self.vsizer.Add(self.lims_hsizer,1,wx.ALIGN_CENTER_HORIZONTAL)
        
        
        self.fit_order_sizer = wx.BoxSizer(wx.HORIZONTAL)
        self.fit_order_sizer.AddSpacer(5)
        self.fit_order_txtbox = wx.TextCtrl(self, wx.ID_ANY)
        
        self.fit_order_sizer.Add(wx.StaticText(self,wx.ID_ANY, "Fit Order:"),0, wx.ALIGN_CENTER_VERTICAL)
        self.fit_order_sizer.Add(self.fit_order_txtbox,0, wx.ALIGN_CENTER_VERTICAL)
        
        self.vsizer.Add(self.fit_order_sizer,1,wx.ALIGN_LEFT)
        
        self.apply_button = wx.Button(self, wx.ID_ANY, "Apply")
        self.vsizer.Add(self.apply_button,0,wx.ALIGN_BOTTOM|wx.ALIGN_RIGHT)
        
        
        wx.EVT_BUTTON(self, self.get_lims_button.GetId(), self.on_get_lims)
        wx.EVT_BUTTON(self, self.apply_button.GetId(), self.on_apply)
        wx.EVT_CHECKBOX(self, self.autofit_check_box.GetId(), self.on_autofit)
        self.SetSizer(self.vsizer)
        self.vsizer.Fit(self)
        self.SetAutoLayout(1)
        self.update(self.plot_manager.get_current_scan())
        self.plot_manager.register_callback(self.update)
        
    def on_apply(self,evt):
        wx.BeginBusyCursor()
        scan = self.plot_manager.get_current_scan()
        old_key = scan.ica_time
        scan.manual_bkgd = True
        scan.bkgd_lowlim = float(self.lower_lim_txtbox.GetValue())
        scan.bkgd_highlim = float(self.upper_lim_txtbox.GetValue())
        scan.bkgd_fit_order = int(self.fit_order_txtbox.GetValue())
        scan.calculate()
        self.plot_manager.update(scan, old_key = old_key)
        self.plot_manager.scan_selector.update_plot()
        wx.EndBusyCursor()
    
    def update(self,scan):        
        
        self.fit_order_txtbox.ChangeValue(str(scan.bkgd_fit_order))
        if not scan.manual_bkgd:
            self.autofit_check_box.SetValue(True)
            self.on_get_lims(None)
            self.set_non_editable()
        else:
            self.autofit_check_box.SetValue(False)
            self.lower_lim_txtbox.ChangeValue(str(round(scan.bkgd_lowlim,2)))
            self.upper_lim_txtbox.ChangeValue(str(round(scan.bkgd_highlim,2)))
            self.set_editable()
    
    def on_autofit(self,evt):
        if self.autofit_check_box.IsChecked():
            self.set_non_editable()
            wx.BeginBusyCursor()
            scan = self.plot_manager.get_current_scan()
            old_key = scan.ica_time
            scan.manual_bkgd = False
            scan.calculate()
            self.plot_manager.update(scan, old_key = old_key)
            self.plot_manager.scan_selector.update_plot()
            wx.EndBusyCursor()
        else:
            self.set_editable()
    
    def set_non_editable(self):
        self.lower_lim_txtbox.Disable()
        self.upper_lim_txtbox.Disable()
        self.fit_order_txtbox.Disable()
        self.get_lims_button.Disable()
        self.apply_button.Disable()
    
    def set_editable(self):
        self.lower_lim_txtbox.Enable()
        self.upper_lim_txtbox.Enable()
        self.fit_order_txtbox.Enable()
        self.get_lims_button.Enable()
        self.apply_button.Enable()
        
    def on_get_lims(self, evt):
        
        lower_lim, upper_lim = self.plot_manager.bkgd_selector.get_lims()
        
        self.lower_lim_txtbox.ChangeValue(str(round(lower_lim,2)))
        self.upper_lim_txtbox.ChangeValue(str(round(upper_lim,2)))


class FileChooser(wx.Frame):
    def __init__(self):
        wx.Frame.__init__(self, None, wx.ID_ANY)
        self.top_panel = wx.Panel(self, wx.ID_ANY)
        self.main_sizer = wx.BoxSizer(wx.VERTICAL)
        
        self.filename_box = wx.TextCtrl(self.top_panel, wx.ID_ANY)
        self.browse_button = wx.Button(self.top_panel, wx.ID_ANY, "Browse")
        self.file_hsizer = wx.BoxSizer(wx.HORIZONTAL)
        
        self.file_hsizer.Add(self.filename_box,1,wx.ALIGN_CENTER_VERTICAL)
        self.file_hsizer.AddSpacer(5)
        self.file_hsizer.Add(self.browse_button,0,wx.ALIGN_RIGHT)   
        self.main_sizer.Add(self.file_hsizer,0,wx.ALIGN_CENTER_HORIZONTAL|wx.EXPAND)
        
        
        self.hdist_txtbox = wx.TextCtrl(self.top_panel, wx.ID_ANY)
        self.hdist_sizer = wx.BoxSizer(wx.HORIZONTAL)
        self.hdist_sizer.Add(wx.StaticText(self.top_panel,wx.ID_ANY,"Horizontal Distance:"))
        self.hdist_sizer.Add(self.hdist_txtbox)
        self.main_sizer.Add(self.hdist_sizer)
        
 
        self.vangle_txtbox = wx.TextCtrl(self.top_panel, wx.ID_ANY)
        self.vangle_sizer = wx.BoxSizer(wx.HORIZONTAL)
        self.vangle_sizer.Add(wx.StaticText(self.top_panel,wx.ID_ANY,"Vertical Angle:"))
        self.vangle_sizer.Add(self.vangle_txtbox)
        self.main_sizer.Add(self.vangle_sizer)     
        
        self.main_sizer.AddStretchSpacer()
        self.ok_button = wx.Button(self.top_panel, wx.ID_ANY, "Ok")
        self.main_sizer.Add(self.ok_button,0,wx.ALIGN_RIGHT|wx.ALIGN_BOTTOM)  
        self.top_panel.SetSizer(self.main_sizer)
        self.main_sizer.Fit(self)
        self.top_panel.SetAutoLayout(1)
        wx.EVT_BUTTON(self, self.browse_button.GetId(), self.on_browse)
        wx.EVT_BUTTON(self, self.ok_button.GetId(), self.on_ok)
        
    def on_browse(self, evnt):
        dialog = wx.FileDialog(self, "Choose input file")
        
        if dialog.ShowModal() != wx.CANCEL:
            self.filename_box.SetValue(dialog.GetPath())
    
    def on_ok(self, evnt):
        filename = self.filename_box.GetValue()
        
        self.Show(False)
        wx.YieldIfNeeded()
        if filename.endswith(".scans"):
            with open(filename, "rb") as ifp:
                scans = cPickle.load(ifp)
            t1 = datetime.datetime.now()
            cProfile.runctx('estimate_ica_error(scans[0],1000)',globals(),locals(),"profile_file")
            t2 = datetime.datetime.now()
            print "time taken = ", t2-t1
            sys.exit()
        else:
            h_dist = float(self.hdist_txtbox.GetValue())
            v_angle = float(self.vangle_txtbox.GetValue())
            conv = 1.0
        
            #run_scan = lambda x: x.calculate()
            
            scans = []
            for i in processing.fullscans(filename, conv=conv):
                #scans need to have at least 4 points in them - otherwise we can't fit a gaussian to them
                if len(i["angle"]) < 4:
                    continue
                scans.append(ProcesssedScan(i["angle"].copy(), i["times"].copy(),i["so2_col_amount"].copy(), i["angle_error"].copy(), i["so2_col_amount_error"].copy(), h_dist, v_angle))
            
            #do the fittings/integrations
            scans = parallel_process(calc_scan, scans)
        
        self.p = PlotManager(scans)    
        ctrl_win = ControlWindow(self.p, os.path.basename(filename)+" - ICA Inspect")
        self.p.show()
        self.Destroy()
        
def calc_scan(x):
    return x.calculate()
        
if __name__ == '__main__':
    
    if len(sys.argv)>=2 and sys.argv[1] == "test":
        import doctest
        doctest.testmod()
        sys.exit()
    
    app = wx.PySimpleApp()

    f = FileChooser()
    f.CenterOnScreen()
    f.Show()
    app.MainLoop()
#    f.ShowModal()
#    
#    filename = wx.FileSelector("Choose file to process.")
#    
#    if filename.endswith(".scans"):
#        with open(filename, "rb") as ifp:
#            scans = cPickle.load(ifp)
#    else:
#        h_dist = 1000.0
#        v_angle = 30.0
#        conv = 1e-4
#    
#        run_scan = lambda x: x.calculate()
#        
#        scans = []
#        for i in processing.fullscans(filename, conv=1.0):
#            scans.append(ProcesssedScan(i["angle"].copy(), i["times"].copy(),i["so2_col_amount"].copy(), i["angle_error"].copy(), i["so2_col_amount_error"].copy(), h_dist, v_angle))
#        
#        scans = parallel_process(run_scan, scans)
#    
#    p = PlotManager(scans)    
#    ctrl_win = ControlWindow(p)
#    p.show()
